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Abstract 

In this article, we show that the project ion-free, snapshot-based, balanced truncation method can 
be applied directly to unstable systems. We prove that even for unstable systems, the unmodified 
balanced proper orthogonal decomposition algorithm theoretically yields a converged transformation 
that balances the Gramians (including the unstable subspace). We then apply the method to a 
spatially developing unstable system and show that it results in reduced-order models of similar 
quality to the ones obtained with existing methods. Due to the unbounded growth of unstable 
modes, a practical restriction on the final impulse response simulation time appears, which can be 
adjusted depending on the desired order of the reduced-order model. Recommendations are given 
to further reduce the cost of the method if the system is large and to improve the performance of 
the method if it does not yield acceptable results in its unmodified form. Finally, the method is 
applied to the linearized flow around a cylinder at Re = 100 to show that it actually is able to 
accurately reproduce impulse responses for more realistic unstable large-scale systems in practice. 
The well-established approximate balanced truncation numerical framework can therefore be safely 
applied to unstable systems without any modifications. Additionally, balanced reduced-order models 
can readily be obtained even for large systems, where the computational cost of existing methods is 
prohibitive. 


1 Introduction 

Many linear dynamical systems, such as the linearized Navier-Stokes equations, are composed of a large 
number of states 0(10^ — 10®), but their behavior is dominated by a much smaller number of modes 
0(1 — 100). Obtaining a low-order model that only retains the dominant features of the system’s behavior 
is of great value in order to understand and modify its dynamics. Specifically in a feedback control setting, 
most controller design methods are only tractable if a reduced-order model (ROM) is available. 

In fluid mechanics, many successful ROMs of the flow field dynamics have been obtained by projecting 
the system equations onto a low-dimensional subspace, composed of a set of particularly relevant modes, 
such as global modes, proper orthogonal modes (ROMs) or balanced modes. Global modes are ranked by 
their damping rate, so projecting the dynamics onto the least stable modes is a natural choice. Akervik 
et al. [3] and Henningson et al. [23] successfully applied this strategy to damp global oscillations in a 
shallow cavity using an LQG controller. However, it is common for the dynamics to be dominated by the 
non-normal interaction between (potentially highly damped) global modes, in which case a large number 
of modes may be required to obtain an acceptable ROM. A more strategic selection of global modes 
can result in better performance in some cases [laiia but finding a set of robust selection criteria is 
not straightforward m- Another weakness of this method is that identifying a large number of (highly 
damped) global modes may be prohibitively expensive. 

Alternatively, proper orthogonal modes provide an optimally low rank approximation of the state over 
a chosen set of snapshots as they are ranked by their energy content. ROMs of many flow-fields have been 
developed by projecting the dynamics onto these modes. For instance, models of boundary layers [51135], 
channel flows [37], backward-facing steps [551 HHl [IH] , bluff body flows and wakes [3^ [351 [351 Sll SSI SI] , 
forward facing steps [41], and cavities [3T] have all been studied. POD-based ROMs are attractive both 
for the simplicity of the (snapshot-based) method used to obtain them [46], and the intuitive projection 
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basis they provide (which allows retaining most of the energy of the snapshot ensemble). However, there 
is not a strong theoretical justification for using this orthogonal basis, as it depends on the specific set of 
snapshots used to construct it. It is therefore unsurprising that a large number of modes are sometimes 
required to accurately represent the flow dynamics m- Several methods have been developed to improve 
the robustness and reliability of the models, such as regularly adapting the set of snapshots used to form 
the basis [40l [13], or adding a so called shift-mode [35] which acts as a mean flow-correction mode and 
allows transients to be modeled more accurately. 

In this article, we focus instead on ROMs based on balanced modes. These are ranked by their 
dynamical significance to the input-output relationship of the system, and are therefore better suited for 
feedback control purposes than ROMs or global modes by design. ROMs based on balanced modes have 
attracted a significant amount of attention in recent years, and have been used in many flow configurations 
such as the response of the vertical force of an airfoil to a plunging motion [^ , a channel flow [33] , a cavity 
flow E], a boundary layer flow 130, the flow over a flat plate at large incidence [T], over a backward 
facing step |18| . in a three dimensional boundary layer [44) . and over a cylinder [50] . Furthermore, several 
studies have compared the performance of ROMs based on global, proper orthogonal and balanced modes 
[52l im [9] [18] and have concluded that balanced ROMs typically require a much smaller number of modes 
for a given degree of accuracy. 

As the classical (exact) balanced truncation method is intractable for large systems, Moore [32) . 
Willcox and Peraire [S2] , and Rowley [33] , developed an approximate snapshot-based method sometimes 
called balanced POD (or BPOD, see Section [H]) to reduce the computational cost of the balancing 
procedure. When using this snapshot method, approximate balanced modes can be obtained from as few 
as two impulse response simulations (for single-input-single-output or SISO systems), and the snapshots 
do not need to be updated or complemented by a shift-mode. 

A further issue with the classical and snapshot methods is that they were designed for stable systems. 
While a number of studies [13 [3011101 [S3 H] have developed techniques to balance unstable systems (see 
13:31) . only a recent extension of the snapshot method mm has allowed the computation of ROMs for large 
unstable systems. The method was applied to the flow over a cavity m, a flat plate at large incidence 
[I], and a cylinder [SO]. This extension requires the computation of the system’s unstable global modes 
and the projection of the system onto its stable subspace. For large systems, such as three-dimensional 
flows, the cost of this procedure may still be excessively large. 

The aim of this paper is to show that this expensive projection step is in fact unnecessary and 
that applying the unmodified snapshot-based balanced truncation method (designed for stable systems) 
directly to unstable systems yields a converged transformation, which balances the system and can result 
in ROMs of the same quality as those obtained when the projection method is used. Many large- 
scale unstable systems (e.g. three-dimensional flows) for which the projection step makes approximate 
balanced truncation intractable can therefore readily obtain balanced reduced-order models by using the 
projection-free method. 

In Sec. [21 the theoretical framework is introduced and the existing balanced truncation methods are 
outlined. The fact that the projection-free snapshot method can be used for unstable systems is proven 
in Sec. [3] and it is then applied to a representative one-dimensional model system in Sec. [3] where it 
is compared to existing methods. In Sec. [S] the method is applied to a two-dimensional Navier-Stokes 
simulation of the flow over a cylinder to show that it also performs well in this more realistic and 
computationally demanding setup. 


2 Background 

2.1 Balanced truncation of stable systems 

In this section, we introduce the main concepts and methods that are relevant to balanced truncation, 
in particular for stable systems. More details about these methods can be found in standard control 
textbooks (e.g. [S3] or 0)- The standard continuous time, linear time-independent state-space system 
is: 


\ X = Ax Bu, 

\y = Cx, 

where x G are the states of the system, u G C"" are the inputs, y G C"® are the outputs and 
A G B G and C G are three time independent matrices. The dynamics of the 
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system are then governed by the transfer matrix: G{s) = C {si — A) ^ B, where Y{s) = G{s)U{s), s 
is the Laplace variable and U{s) and ^(s) are the Laplace transforms of the input and output signals, 
respectively. 

The balanced truncation approach is based on an analysis of the controllability and observability of 
the system. The controllability of a state is related to the minimum input energy required to reach it 
from a;(0) = 0. The observability of a state Xq is related to the energy of the output signal generated by 
the system starting from x(0) = xq, without any input: u{t) = 0. A state that has a large impact on the 
input-output behavior of the system - i.e. on the transfer matrix - is said to be dynamically significant. 

If a state is unobservable or nearly unobservable, then even if only a small amount of input energy 
is required to reach it (i.e. if it is highly controllable), it will not have a large impact on the output 
signal. Conversely, if a state is uncontrollable or nearly uncontrollable, then a large (or infinite) amount 
of input energy is required to reach it, so only a comparatively negligible part the output energy can be 
due to that state (as long as a;(0) = 0). Balanced truncation therefore aims to create a ROM by only 
retaining the states whose observability and controllability is high, as these are the ones which have the 
largest impact on the system dynamics. In order to identify which states are the most controllable and 
observable, let us define the controllability Gramian: 


H^c(foo) 





and the observability Gramian: 


W'o(too) 





( 2 ) 

(3) 


where ^ is the complex conjugate transpose. These are defined for stable and unstable systems for 
0 < too < + 00 . As too —>■ +00 however, the Gramians converge to constant matrices for stable systems 
but become unbounded for unstable systems. When the limits exist, the following notation is used: 
hm^^_^_i_Qo ITo(too) — bLo and lim^^_^_|_oo bbc(too) — bLc- 

A stable system is said to be balanced when Wc = Wo = where E = diag ([ o'! ... tT„^ ]) and 

E G is a diagonal matrix, where ui > a 2 > ■ ■ ■ > Unx are the Hankel singular values (HSVs) of 

the system. The (non-zero) HSVs are unique and provide an indication of the corresponding state’s 
dynamical significance. In other words, the states of a balanced system are ranked according to their 
joint controllability and observability. 

In general however, systems are not balanced and given a transfer matrix G(s), the realization 
(A, B, C) or internal coordinate system used to define the states x is not unique. In order to balance a 
system it is therefore necessary to find the coordinate transformation T = S~^ G ensures: 

SWcS^ = T^WoT = E^. The transformed (balanced) system is then: 


X = SATx + SBu^ 
y = GTx. 


In order to obtain T and S for stable systems, the converged Gramians are evaluated by solving the 
following Lyapunov equations: 


A^Wo + WoA + G^G = 0 , 

AWe + WoA^ + HRt = 0. 


The balancing transformations are then found using Eq. (15^ . which can be computed from the singular 
value decompositions (SVDs) in Eq. (IHal) and (IHbll : 


Wc = XX^, Wo = ZZ^, 


Z^X = GEVt = 
T = AVlEC^''^ 


'El O' 



0 0 


[ V2^ j 


s = sf^''^ulzK 


(6a) 

(6b) 

(6c) 
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Note that X and Z are not unique. The balanced system can then be decomposed as follows: 


Xi 


1 

Ai2 


Xi 


’ Bi ' 


= 





-1- 


X2 


_1 

A 22 


X2 


. ^2 . 


u, 







Xi 

y = 

Gi G 2 






X2 


and as the transformed states are now ranked by dynamical significance, the less observable and control¬ 
lable states X 2 can be truncated to give the following ROM: 


{ Xi = AiiXi -I- Biu, 

y = Cixi, 

which is also a stable balanced system. An attractive property of balanced truncation is that an a priori 
upper-bound on the ROM error exists: ||G — Gr.||oo < where r is the number of states that 

have been retained in the ROM, H-Hoo is the 'Hoo norm, and Gr{s) = Gi(s/ — is the reduced 

transfer matrix. 

Unfortunately, if the system matrices are not explicitly available or if the system dimension is too 
large, this technique cannot be applied directly. Moore [32], Willcox and Peraire [52] and Rowley 
therefore introduced an approximate balanced truncation method sometimes referred to as balanced POD 
(BPOD), which is based on snapshots from the impulse response of the primal system (|T]) and the adjoint 
system 0, also assumed to be available: 

fi = Atz + Gt^;, 

\w = B^z. 


In this approach, snapshots from the primal impulse responses x{tck) = Xk = and the adjoint 

impulse response z{tok) = Zk = ^re stacked into matrices in order to find an approximation 

for the Gramians in the form Eq. (|6ap . by defining X and Z in the following way: 

X=[xiy/^i ... XN^^/^^]^W,^XX\ (8a) 

z = [ Zi . . . ZNo \/ ^oNo ] ^ Wo ^ZZ\ (8b) 

where there are primal snapshots and No adjoint snapshots taken at discrete times tck and tok 
respectively. Sck & R and 6ok & R are the associated quadrature coefficients corresponding to a chosen 
numerical integration scheme. As a result, XX"^ and ZZ‘^ are discrete versions of the continuous definition 
of the Gramians in Eq. Q and Eq. dSj), respectively. Note that Xk S and Zk € and that 

the snapshots are not necessarily equally spaced in time. Forming X and Z in this manner therefore 
allows T and S to be computed using Eq. (1^ . but at a reduced cost since it avoids finding the solution 
of the Lyapunov equations m and also only requires finding the SVD of a {NcUu x NoUy) matrix instead 
of the three {rix x rix) matrices in Eq. (I6al) and This can represent significant savings if NcUu rix 
and NoUy <C rix, which is usually the case in computational fluid dynamics. Additionally, “squaring up” 
the matrices X and Z to form the Gramians is detrimental to the accuracy of the results [32] . 


2.2 Balanced truncation of unstable systems 

In this section, existing methods related to the balanced truncation of unstable systems are discussed. 
For stable systems it is straightforward to show that the balancing transformation T converges and that 
it diagonalizes and equalizes the Gramians as too -foo, since the primal and adjoint states go to zero. 
If the system is unstable however the state becomes unbounded for large too- ITc(too), Wo (too) as defined 
in Eq. @ and dSj) respectively are therefore also unbounded and Wc and Wo are not defined. 

Nevertheless, Ghiu Kenney and Hewer Therapos [15], and Al-Saggaf H] showed that despite 
Eq. ([5]) and dSj) being unbounded, Eq. dSj) still have solutions if and only if Xi + Xj ^ 0, where Xi and Xj 
are any two eigenvalues of A. When this is the case, a balancing transformation can be obtained if WcWo 
is similar to a real diagonal matrix. However, this is problematic when Eq. m does not have a solution 
and when the balancing transformation does not exist. Furthermore, the resulting reduced-order models 
are not always of satisfactory quality [51] . 
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An alternative method developed by Meyer [30] is based on the co-prime factorization of the trans¬ 
fer matrix. Nett [33] showed that G{s) = N{s)M{s)~^, where N{s) = C{sl — A)~^B and M(s) = 
/ -I- K{sl — A)~^B is a right co-prime factorization of G{s) if A = A + BK is stable. Meyer and 
Franklin m showed that ii K = —B^P, where P is the solution to the algebraic Riccati equation; 
PA + Alp — PBB^P + C^C = 0 then the co-prime factorization is normalized. A realization of the sta¬ 
ble system [ N\s) M^s) ] ^ is then given by B, ( i balanced and 

truncated as in Sec. 12.11 The stabilizing feedback K and output augmentation can then be undone to 
retrieve the ROM. An issue with this method is that it does not simplify to the standard method when 
the plant is stable. 

Zhou |S3] developed a related method by using the frequency-domain definition of the Gramians: 

1 -1 
W^f = — (fw/-A)"^PPl (-fwl-Al) dw, 

J — OO 

1 r+oo ^ 

Wof = — (-iwi-Al) C'<C{iojI-A)~'^duj, 

J — OO 

where Wcf and Wof are now well-defined as long as A has no eigenvalues on the imaginary axis. In 
the stable case, Parseval’s theorem can be applied to show that Wcf = Wc and Wof = Wo- It is then 
possible to evaluate Wcf and Wof if a transformation that decouples the stable and antistable dynamics 
of the system is available. However, they are usually instead evaluated by hnding the stabilizing so¬ 
lutions to (A -f BK)Wcf + Wcf {A + BK^ + PPt = 0 and Wof {A + LG) -b (A -f LC)'<Wof + C^C = 0, 
where K = —B^P and L = —YG. Here P is the stabilizing solution to the algebraic Riccati equation 
PA + Alp — PBB^P = 0 and Y is the stabilizing solution to the algebraic Riccati equation AY -|-FAl — 
YC^CY = 0. The balancing and truncation procedure from Sec. I2.1l is applied once these “generalized” 
Gramians are known. 

A snapshot-based extension to this method was developed by Dergham et al. [18] , who applied it to a 
rounded backward facing step and cavity flow. Here, the snapshots are defined in the frequency-domain 
rather than the time-domain: the primal system snapshots are of the form X (w) = {iujl — A) ^ B and 
of an analogous form for the adjoint snapshots. Each snapshot therefore has to be evaluated by inverting 
a large matrix explicitly, but all the snapshots can potentially be evaluated simultaneously in parallel. 

In the methods described above, the Hankel singular values corresponding to unstable modes are not 
necessarily larger than the ones corresponding to stable modes, so there is no guarantee that unstable 
modes will not be truncated. This is often not a desired property for control design purposes. If on 
the other hand the transformation that uncouples the stable and antistable dynamics is available, such 
that G{s) = Ga{s) + Gs{s), then an alternative approach, for instance suggested by Enns [5D], is to only 
balance and truncate the stable part of the system Gs{s) and simply add the unbalanced antistable 
dynamics back in such that the ROM transfer matrix is: Gr = Ga{s) +Gs{s). where " refers to the 
balancing and truncation procedure, thus conserving all the unstable modes. 

Unfortunately, if only a time-stepping code is available and if the system dimension is large, none of 
the methods described above can be used. In order to tackle this issue, an extension to the snapshot 
method introduced in Sec. l2.1l was developed by Barbagallo et al. [TT] and Ahuja et al. [Tj. This extension 
is closely related to the approach proposed by Enns [20] as it is based on the separation of the stable and 
antistable dynamics of the system. Although it may be expensive to compute the full stable-antistable 
uncoupling transformation, it may still be possible to use an Arnold! procedure to identify the right and 
left antistable eigenspaces Pa and Qa respectively (scaled such that Ql^Pa = I). The primal and adjoint 
systems can then be projected onto their respective stable subspaces, by defining the projection matrix 
V = I-PaQl 


is = PAVxs + PBu, 

Vs = GVxs, 


and 


\zs = + ptc-ti;, 

= B^^V^Zs. 

This projected system is then balanced and truncated using snapshots, as described in Sec. 12.11 The 
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antistable dynamics are finally added back in, and the final reduced-order model is: 


Xr = ArXr “I" BrU, 

y = CrXr, 


where 




■ QlAP, 

0 ■ 

^ = 

■ QlB - 

0 

A 


B 


[CPa C], 


where, as above • refers to the projected system’s balanced and truncated matrices. Like in the method 
proposed by Enns [5D] , this procedure does not balance the unstable subspace (unlike the ones developed 
by Kenney and Hewer Meyer [30], and Zhou [M]), but guarantees that unstable modes are not 
truncated. For fluid systems, the dimension of the unstable subspace is often 0(1 — 10) so this method 
often only requires computing a few eigenvalues and eigenmodes. However an Arnold! package is not 
always available and for large systems (e.g. three-dimensional flows), finding even a few eigenmodes may 
still not be tractable. A further issue with this method is that if the system is poorly conditioned, it may 
be difficult to evaluate the projection matrix itself. 

In this article we propose an alternative method that is projection-free and does not require evaluating 
any global modes. As briefly mentioned in |^, one can choose to simply use the standard method for stable 
systems, but based on finite-time Gramians, which are expected to approximate the impulse response 
over the chosen time interval well-enough. This article aims to show that if this finite time interval is 
long enough, the balancing transformations actually converge and hence balance the Gramians for any 
further time. These transformations can therefore be considered to be true balancing transformations, 
often obtained at a fraction of the cost of the methods described above. As a result, we show that the 
projection-free, snapshot-based balanced truncation method is directly applicable to unstable systems. 


3 Projection-free balanced truncation of unstable systems 

3.1 Theoretical justification 

The goal of this section is to show that the projection-free snapshot-based balanced truncation method 
(BPOD) introduced in Sec. l2.1l can be used even in the unstable case, despite the fact that the Gramians 
are unbounded. In order to complete the proof we set out to prove the following statements: 

1. The transformations T and S as defined in Eq. (13^ converge to constant matrices, even for unstable 
systems. 

2. The controllability and observability Gramians are balanced by the converged transformations for 
any sufficiently large too, despite not converging to constant matrices. 

These two statements are proven in Proposition 2 and Proposition 3 respectively and ensure that the 
converged T and S matrices can be considered to be balancing transformations for the unstable system. 
These proofs are based on Proposition 1, which states that as too +oo, the unstable singular vectors 
and values of the matrices X and Z can be identified explicitly. In particular, all the left singular vectors 
tend to constant vectors. The proofs outlined here correspond to the case where all eigenvalues have 
distinct real parts. Appendix [^ deals with configurations where some of the modes have the same growth 
rate. 

Let (A, B, C) be a minimal realization, where A = PAQ\ and P^Q = /, where P = [ pi ... Pn^ ] , 
A = diag ([ Ai ... A„„, ]), Q = [ ... ], and Re(Ai) < ... < Re(A„„,). The primal and ad¬ 

joint impulse responses are defined as in Sec. 12.11 and can respectively be written: 

x{tck) = Xk = Pe^^*'=^''Q^B, 
zitok) = = 

The approximation of X and Z from the impulse response snapshots defined in Eq. m can be written: 

X = P/3t = Z = = Uo^oV^, 
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where /3 = [ ... /3n^ ] and C = [ ] and: 


A = 


e^ = 




where the superscript * refers to the complex conjugate. 


Proposition 1: For large too, the singular values and vectors of X (Z) corresponding to 
unstable singular values in Ec (Eq) converge to: 


{ ^ci — ^ci—iPi II ^C2—iPz II : 

aci = ||rcj_ipj||||^j||, 

Vci = /3i||Air\ 

^oi — ^oi—1 II ^oz—iPz II : 

ao^ = ||To,_igi||||^i||, 

Voi = ■fz|l?iir\ 

respectively, where H H is the Euclidean norm of a vector and: 


Pci — i- [ ^cl ■ . • ^cz 


Toi — f [ ^ol . ■ ■ ^oi 


r „t 1 


1 


(9a) 


(9b) 


(10a) 


(10b) 


Proof. We will use a proof by induction and choose the induction hypothesis (I„) to be the fact that the 
proof holds for alH < n. 

For the base case (Ii), if Re(Ai) > 0 then the first mode pi is unstable and: 

lim Xk = 
fc—>■ + 00 

lim X = pi+ = zicio'ciulti, 

too—^+00 

since this is just rank 1. iXx) therefore holds, up to a unit norm factor 

{uci ^pilbill'^e*®^!, 

<^cl = ||pi|||bl||, 

[vci = ^i|bi||“b-*®“i, 

where 9ci is a real scalar. The first singular vectors and value are thus uniquely defined, up to 
which we can always choose to be equal to 1, so we will make this assumption without loss of generality 
for the remainder of this paper to simplify the notation. As a result, the direction of the first left singular 
vector Uci converges to that of the first eigenvector of the system pi (which is constant). 

In the inductive step, we prove that if (I„) holds for some rank n, then (In+i) also holds. The 
transformation matrices Td defined in Eq. (I10a|l project out all the left singular vectors Ucj for all j < i, 
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i.G. TQ^Ucj — 0, TQjUf^i — Uci- 


— [ '^cn+l ■ ■ • '^cux ] 


^cn+1 


^cn+1 


— [ T^cnPl • ■ • '^cnPn ] 


Pi 


Pi 


[ '^cnPn+l ■ • ■ T^cnPrix ] 


*" /^l+i 


/3L 


Thus, if Pn+i is an unstable mode and recalling that we assumed that Re(Ai) < ... < Re(A„^), we have: 


lim. TcfiX — [ TqyiPi ■ • ■ TcnPcn 1 

too—^ + 00 '' 


?! ^ 


Pi 


H“^cnPn+l/5yj_|_i ■ 


Now for 0 < i < n, Ten can be written: 



/ 

’ ' 

\ 

T — 

^ cn — 

I 1^ Uci ■ ■ • Ucn ] 




1 

Tit 

L ^cn J 



T 

C.l. — 


CZ —1 5 


and since we are assuming {Xn) holds for rank n\ 



/ 

' X 

\ 

TenPi — 

I [ '^cz ■ ■ • '^cn ] 




1 

L ^cn J 

) 


ll^CZ—II - O 5 


lim Tq^X — TqyiP^-^iP — '^cn+l^cn+l^cn+l ’ 

too—^ + 00 

and hence {Tn-\-i) holds too. This completes the inductive step and along with the base step concludes 
the proof by induction. The singular values and vectors corresponding to unstable modes are therefore 
given by: 

^ci — Tci—1^2 II ?ci_ iPi II , 
ac^ = ||Tci_iPi||||/3i||, 

^Vci = Pi\\Pi\\~^. 

The singular vector Ud is therefore pointing in the direction of the component of pi that is orthogonal to 
the subspace defined by [ pi ... Pi-i ]. The procedure that identifies the left unstable controllability 
singular vectors is therefore essentially a Gram-Schmidt process. An analogous derivation starting from 
the adjoint impulse response leads to: 

^oi — II II , 

aoi = ||Toj_ig*||||^i||, 

v„^ 

completing the proof of Proposition 1. □ 

It is clear that the left singular vectors and singular values of X and Z, when restricted to the stable 
subspace also converge to constants as toe —t + 00 , (this is the basis for the snapshot method in [TT] and 
DP)- Therefore all stable and unstable left singular vectors of X and Z tend to constants, i.e. Uc and Uq 
converge to constant matrices. 

Proposition 2: The balancing transformations T and S converge to constant matrices for 
large too- 
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Proof. Let us define the Hankel matrix Mh, noting that the SVD of the Gramians can be approximated 
with snapshots for any finite too: hLc(too) ~ XX\ Wo(too) ~ ZZ"^: 

Mh = Z^X = = KSof/tf/oSoy;. (11) 


Using a similar reasoning to the one in the proof of Proposition 1, the (i + l)th set of Hankel singular 
values and vectors corresponding to each unstable mode of the system can be identified by projecting Mh 
onto the subspace that is orthogonal to the left and right singular vectors corresponding to all unstable 
modes j such that j < i. By using a similar proof by induction to that in Proposition 1 (see Appendix [El 
for the full derivation) and by defining the following two transformations: 


Tu = I - \ ui ... Ui 


Tri — L [ Vi ... Vi 


r 1 


V 


it can be shown that as too —>■ +oo: 




Vi = Vro.. 


(12a) 


(12b) 


(13) 


Because the singular values corresponding to unstable modes tend to infinity for large too, we can separate 
the stable and antistable parts of X and Mh (denoted with the subscripts s and a respectively) as follows: 


lim X 

too —>- + 0O 


lim Mh 

too —>-+00 


Now, Eq. (fTEI) implies that = /, 

T, as defined in Eq. (IHEl becomes: 


lim T = 

too ^+cx:) 


Eca 0 

r ut 1 

* ca 

0 Ecs 

ut 

cs 


= [ Uoa Uos ] 

= [Ua Us] 

VfoYs = 0, and VfsVa = 0. As a result, the transformation matrix 


\ E„ 

0 

r ut 1 

a 

O 

_1 

Eg _ 

-1 

S”” 

_1 


Uoa^oa^a Uos^osVlVs^s 


[Ta Ts]. 


(14) 


The tth column of the converged Ta is therefore Tai = Ud The ratio of the singular values can be 
shown to tend to a constant (this fact is proven in Appendix [C]) so the matrix T converges to a constant 
matrix for large too (since Tg must also converge). An analogous argument can be made to show that 
S converges to a constant matrix too, where lim(^_>+oo <5 = [ S'^ 'S'l ] , where Sa = Ea ^^‘^Z^oaUlg^ and 
Ss = This completes the proof of Proposition 2. □ 


Proposition 3: The converged balancing transformation T {S) balances the Gramian Wo 
(Wc) for sufficiently large too- 

Proof. If a (converged) transformation T{ti) is found, corresponding to a given set to snapshots with 
the final snapshot taken at too = ti, we would like to check that it diagonalizes and equalizes the Grami¬ 
ans Wc{t 2 ) and Wo(t 2 ), corresponding to a different (but also sufficiently large) set of snapshots such 
that too=t 2 - Using the notation =Z(t2)^-A(ti) = {7(2i)S(2i)U^2i)> transformed observability 

Gramian becomes: 


T^WoT = r(ti)tlTo(t2)T(ti), 

« (E-;f U(\,)Xt(U)) Z{h)Z{t2)^ (x(U)U(n)E-/,f) , 
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and V( 2 i) can be obtained from the matrix Mh{ 2 i) = ^4(^2) [^o{t 2 )U^{t 2 )Uc{ti)Y,c{ti)] By fol¬ 

lowing the same procedure as in the proof of Proposition 2, we obtain: 

Ua{21) = Voa{t2), K(21) = 14a (^l) = 14(11), (16) 

where, as above, the a subscript refers to the antistable part of the matrix for large too (i-C. large ti and 
t 2 ). For the stable part, given Eq. (ITT)1) and the fact that the stable modes decay for large ti and t 2 , any 
additional snapshots cannot modify these subspaces so 14 ( 2 i) = 14(ii)- Finally: 

^ ^ SCA)S^21): (17) 

which is clearly diagonal. An analogous proof can be applied to the controllability Gramian and leads to 
S'(ti)bFc(t2)5'^(ti) ~ ^(ii)^( 2 i) ~ F(ti)lWo(t 2 )F(ti). This completes the proof of Proposition 3. □ 

In Appendix|3 it is shown that the same conclusions hold for unstable and marginally stable repeated 
eigenvalues, as well as marginally stable complex conjugate pairs. For unstable complex conjugate pairs, 
only the two-dimensional subspace spanned by the two corresponding balanced modes converges. This is 
all that we require, as it is usually not desirable to truncate only one of the two modes corresponding to an 
unstable pair of complex conjugate eigenvalues. In this case, the two columns (rows) of the transformation 
T (S) that define this converged two-dimensional subspace oscillate as too is increased. The transformed 
Gramians are therefore equal and diagonal except for one (2 x 2) block along the diagonal for each unstable 
complex conjugate pair. Let us emphasize here that despite the transformation not becoming constant 
and the Gramians not being strictly diagonalized and equal, any T (and S) corresponding to a sufficiently 
large too can be used as the balancing transformation, as long as the 2D subspace corresponding to the 
unstable modes has converged. In this case the unstable complex conjugate mode pair is balanced as a 
whole. 

3.2 Practical considerations 

As illustrated in sections 2] and [SI the projection-free snapshot-based balanced truncation method can 
often be applied directly to large unstable systems, just as in the stable case. Nevertheless, some modifi¬ 
cations to the method that can lead to significant improvements in the quality of the ROMs and/or the 
computational cost of the method are outlined in this section. These may be required for particularly 
challenging systems. 

3.2.1 Final simulation time and sampling intervals 

For both stable and unstable systems, the sampling intervals must be small enough to capture the 
highest frequencies of interest in the flow field and the final simulation time too must be large enough for 
all stable modes to decay. For unstable systems too must also be large enough for the impulse response to 
be dominated by the unstable modes at the end of the simulation, thus allowing Ta and Sa to converge. 

Although there is no theoretical upper-bound on too, in practice as too -boo any initial transients 
eventually become so negligible compared to the long term response that the information related to modes 
that are more stable than the dominating unstable mode(s) may be lost. This can result in an inaccurate 
identification of the corresponding balanced modes. Both the upper and lower bounds on too are clearly 
problem-dependent and the trade-off between these two limits is investigated in further detail in Sec. 21 

3.2.2 Improving the accuracy of the method 

If several unstable modes have nearly identical growth rates or if it is necessary to identify slowly decaying 
modes, the lower bound on too may be higher than its upper bound (see 13.2.11) . It may therefore not 
be possible to obtain sufficiently accurate information about the less unstable parts of the system from 
the impulse response. One way to work around this is to use the fact that for sufficiently large too the 
most unstable mode(s) are still identified accurately. Another set of impulse response simulations can 
then be run with these modes projected out in order to identify the more stable modes. The procedure 
can potentially be repeated until each unstable mode has been identified and projected out. In this 
case the method is clearly not projection-free anymore. It becomes more similar to the one suggested by 
Barbagallo et al. [a and Ahuja et al. [T], although here the entire unstable subspace does not necessarily 
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need to be projected out for the method to work. Additionally, it still avoids the need for an Arnoldi 
solver and theoretically still yields the same transformation T as the projection-free method (and hence 
the unstable subspace is still balanced). 

An alternative approach is inspired from the work of Morgans and Dowling |d3] and Illingworth, 
Morgans and Rowley E3' the unsatisfactory ROM obtained from the impulse response can be used as 
a first approximation to the system (with again the most unstable mode(s) accurately identified). Using 
this ROM, a controller aimed at suppressing the most unstable mode(s) can be designed and the resulting 
closed-loop impulse response can be used to deduce the corresponding open-loop dynamics and improve 
the initial ROM. 

3.2.3 Large systems and Gaussian quadrature 

In order to generate the balanced ROM, snapshots are typically recorded at regular time intervals. The 
matrices X and Z as defined in Eq. (|S]) - which are used to approximate the Gramian integrals in Eq. @ 
and m - are formed by using Newton-Cotes quadrature weights of a selected order (e.g. trapezoidal, 
Simpson or Boole rule). However, the number of snapshots required can be significantly reduced by 
selecting a Gaussian quadrature rule instead (e.g. Gauss-Legendre quadrature), where the snapshots are 
not equally spaced in time. 

Since the final simulation time is usually not known a priori, it may be preferable to use a “composite” 
quadrature rule. In this case each impulse response is divided into several time windows (of potentially 
different lengths). If there are Ni snapshots in the ith time window, the integral of this window can 
then be evaluated independently from the rest of the impulse response, using a (fVi)-point Gaussian 
quadrature. With this piecewise integration method, it becomes straightforward to add an additional 
time window to increase the final simulation time. It is then also possible to optimize the distribution of 
the snapshots across the full simulation, for instance, if a higher-order quadrature is required in a specific 
time window. 

If the system dimension is so large that storing each snapshot is computationally demanding or if 
the number of snapshots required is so large that it results in an excessively expensive SVD (for three- 
dimensional flows for instance), then it may be necessary to use Gaussian quadrature to make the balanced 
POD method tractable. 


4 Application to a spatially developing unstable system 

4.1 Linearized complex Ginzburg-Landau equation simulation setup 

In order to check that the application of the approximate balanced truncation method to unstable systems 
results in satisfactory ROMs for large-scale problems, we first apply the method to the linearized complex 
Ginzburg-Landau equation. This one-dimensional system exhibits spatially developing behavior and 
instabilities, which are representative of typical flow fields. The linearized complex Ginzburg-Landau 
equations (I18al) and the corresponding adjoint equations (IlSbIl are: 



( d 92 





dq+ 


= = 


dx 


dx"^ 


+ ti{x)* q 


(I 8 a) 

(18b) 


where p,{x) = p,o — + pL 2 X^ 12. Note that now x is the spatial variable, while q is the system state. 
The complex parameters v and 7 characterize the convection and diffusion/dispersion in the system 
respectively, while /i is related to the exponential growth of instabilities, /tq can be seen as being similar 
to the Reynolds number in the Navier-Stokes equations, as it is used to determine the nature of the global 
stability of the system. A large value of pL 2 corresponds to a large degree of non-parallelism, while a small 
value of pL 2 and a large value of v result in a strongly non-normal flow. Finally c„ is the most unstable 
wavenumber in the flow-held. 

The results in this section were obtained by using the code developed by Bagheri [S] with a similar set 
of parameters to the supercritical (globally unstable) system considered in [^, as summarized in Table 
[TJ In order to demonstrate the method’s ability to deal with several unstable modes however, ^0 was 
set to 0.57 in order to obtain a two-dimensional antistable subspace. In this code, the forward problem 
uses a spectral Hermite collocation method, where state q G is evaluated at the roots of the rixth 
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Parameter Value 


rix 

220 

Mo 

0.57 

M2 

-0.01 

V 

2 + Q.2i 

7 

l-i 

Xj, Xu 

±10.7 

C’U 

0.2 


Table 1: Ginzburg-Landau equation simulation parameters 


Hermite polynomial Hn^{bx), and the parameter b is chosen to obtain an accurate approximation of the 
continuous problem. The adjoint equations have been obtained from the discretized forward equations 
using a scaling matrix M, so that the energy of the discretized state, defined using the inner product 
q^Mq, approximates the energy of the continuous state. Therefore, the adjoint state-space matrices 
are of the form (A+, i?+, C+j = as opposed to simply (see the 

Appendix of [S] for more details regarding the discretization of the code). 

In the present work, the system was set up as a single-input-single-output (SISO) system: the input 
(actuator) has a narrow Gaussian distribution, centered at xj, the upstream limit of the region where 
instabilities are able to grow (“branch I” in m)- The output (sensor) is defined by the same Gaussian 
function, but centered at a;//, the downstream limit of the growth region (“branch 11”). The discretized 
system has 220 states, corresponding to a spatial extent of [—85,85]. 

The Ginzburg-Landau equation has been often used to model fluid instabilities (reviews on this topic 
such as [241115] frequently demonstrate important concepts with this one-dimensional model), and is 
a common test case for flow-control and model reduction studies of convectively and globally unstable 
flows because its behavior is often representative of the Navier-Stokes equations but at a much lower 
computational cost (e.g. [3S1I1H1IIS])- For more details about the Ginzburg-Landau equation and related 
studies, the reader is referred to [3]. 

4.2 Numerical results 

In this case, the cost of a simulation is low, so there was no need to use a Gaussian quadrature scheme 
suggested in Sec. 13.21 and instead a Boole rule quadrature scheme was chosen. The number of snapshots 
used was chosen such that too/N ps 0.05. The balancing transformations and ROM were then obtained 
as described in Sec. [Q and Sec. [31 As there are only 220 states in the full system, it was also possible 
to compute the transformations and ROMs based on the various procedures described in Sec. 12.21 For 
clarity however, we choose to compare our method to the snapshot-based projection method mm and 
what can be referred to as its “exact” equivalent [20j, where the stable subsystem is balanced exactly. 

4.2.1 Comparison of the system’s behavior with theory 

As the conclusions of Sec.[3]are based on the theoretical limiting behavior of unstable systems, it is crucial 
to investigate the extent to which the key steps of the method yield the predicted results. In this section, 
we therefore compare several matrices obtained from the simulations with their theoretical counterparts. 
When the matrices considered here have a dual analog, which behaves in the same way, we choose to 
focus only on one of them for brevity. 

The first result that is checked here is part of the conclusion of Proposition 1, i.e. that limt^_>_|_oo Ud = 
Uci, where Ud = Td-iPi\\Td-iPi\\~^ for unstable modes (there are two in this case), as stated in Eq. (|9a|l . 
We therefore plot the first 10 columns of jt/JCd, with C/c = [ Uci Uc 2 ], which we expect to be equal 
to [ ]. Figured] shows that indeed, for large enough too, the first unstable singular vector of 

the controllability Gramian Ud points in the direction of pi and the second singular vector Uc 2 points in 
the direction of the component of p 2 that is orthogonal to pi. An analogous conclusion can be drawn for 
Uo and Eq. daS). 

The next result of interest is an intermediate conclusion of Proposition 2: Eq. states that the 
left and right unstable Hankel singular vectors tend to the right unstable singular vectors of Z and X 
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(a) 



2 4 6 8 10 


(b) 

Figure 1: (Color Online) First 10 columns of ICj'C/d for too = 20 (a) and too = 80 (b). The green dashed 
line separates the stable and antistable parts of the matrix. The color scale is from 0 (black) to 1 (white). 



Figure 2: (Color Online) Top left (10 x 10) elements of \V^V\ for too = 20 (a) and too = 80 (b). The green 
dashed lines separate the stable and antistable parts of the matrix. The color scale is from 0 (black) to 
1 (white). 
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Figure 3: (Color Online) Top left (10 x 10) elements of \U^T\ for too = 20 (a) and too = 80 (b). The green 
dashed lines separate the stable and antistable parts of the matrix. The color scale is from 0 (black) to 
1 (white). 



respectively. We therefore expect the following to hold: 


lim V^V = lim 

ioo—^ + 00 too—^ + 


■ KlK ■ 


■ /2 x 2 0 

. K+a _ 


0 KtF. _ 


Indeed for large enough too, Fig. [2] shows that this is approximately true although at too = 80, residuals 
still appear in the cross-terms V^^Vs and V^^Va- 

We now turn our attention to the transformation matrix, where we theoretically expect the direction 
of the unstable balanced modes to tend to that of the unstable left controllability singular vectors Uca- 
In Fig. [21 we plot \UlT\, where Ti are the normalized balanced modes Ti = Ti\\Ti\\ ^ so each column is 
in the same direction as the corresponding column of t/lT: 


lim C/JT = lim 

too^ + OO too—>-1-00 




= lim 

too—>• + 00 




EcaEa 

0 


- 1/2 


We therefore expect: 


U^csTs 




- 1/2 


lim U^T = lim 

too—>- + 00 too—>• + 00 




o 

1—1 

. t/t+o C/t+, _ 


[ 0 UlTg J 


However, as seen in Fig.|21 V^V is not fully converged at too = 80 and hence the only part of the expected 
behavior that clearly appears in Fig. [3] is \U^gTa\ = 0. By recalling that Ucs spans the component of the 
stable subspace that is orthogonal to the antistable subspace, this can be interpreted as the fact that the 
unstable balanced modes span the same subspace as the unstable global modes in practice. However, it 
does not seem that the stable balanced modes are orthogonal to the unstable subspace at too = 80 since 
or that each unstable balanced mode is in the direction of the corresponding controllability 
singular vector since ^ I. These observations can be explained by considering the structure of 

\UIT\ when V^V is still converging and is instead of the form: 


VjV = 


I ^aa 
^sa 


^as 


for some potentially small Caa, Cos, and ega matrices of the correct dimensions: 


[/+ = 


Eca + £as) Eo ^ EcaCasEs ^ 


^-1/2 






(19) 


Given that Scs and Ss converge to constants, while limt^_>+oo E^ = -|-oo and limj^^^+oo ^ca = +oo, it 
is not surprising that limt^_>_|_oo EcsCsaEa ' = 0 and that unless all the elements of Cas are exactly 0, 
limi^_).+oo EcaCasEs ^ 0. Similarly, the off-diagonal elements of Eca (I + Sas) Eo are of the form: 
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Figure 4: Convergence of the balancing transformations T and S. (a) Evolution of the energy of the 
forward impulse response (dashed line) and adjoint impulse response (solid line) states with respect to 
time, (b) Rate of change of the transformation matrices with respect to the final simulation time too used 
to compute them, for the first Ux columns of T and rows of S where At = 0.2. 


_1 /o 

o'cai^aaijO'aj ^ some potentially small scalar Caaij- For large too, if * > j (below the diagonal), then 
Ccai ^ whereas if t < j (above the diagonal), we have Ucai ^ This therefore explains the 

_1 /o 

upper-triangular nature of Uca {I + Cas) Sa ' . 

The behavior of the converging matrices considered up to this point can therefore readily be ex¬ 
plained. We now wish to check that the transformation matrices converge as expected in practice. In 
order to do this, the transformation matrices T and S are computed for a range of too values, and the 
rate of change in the transformations over a fixed time interval At: ||T(too) — T{too — At)\\F/At and 
||<5'(too) — S{toc, — At)\\F/At are shown in Fig.^D (H-Hi;’ is the Frobenius norm of a matrix). Clearly both 
matrices converge to constants for large too, as expected despite the exponential growth of the energy of 
the primal and adjoint states, shown in Fig. |4^. 

It is important to note that although T and S converge, the convergence of the direction of the 
balanced modes also needs to be checked. We can therefore plot 1—|Ti(t)'lTi(t — At)| for each normalized 
balanced mode Ti to check that the direction of the balanced modes also converges. Figure [S] shows that 
even for large too, the most unstable mode converges, while the more stable modes eventually become 
so negligible compared to the most unstable mode that they lose their accuracy. Not surprisingly, the 
less dynamically significant the mode, the lower the value of the maximum too before it diverges. This 
behavior is expected as was mentioned in Sec. 13.21 and the consequences of it for the quality of the ROM 
are further investigated in Sec. 14.2.21 and in the next paragraph. 

Finally, we check that the Gramians are indeed balanced by the transformation. By definition, 
Gramians computed using a given set of snapshots are exactly balanced by transformations T and S, which 
were computed using the same set of snapshots. However, we wish to investigate how well the “converged” 
transformation balances Gramians computed using another set of snapshots. In order to do this, recall 

from Eq. (fT^ that the transformed observability Gramian is (^(ii)^(2i)) ^(21) (^(ii)^(2i)) 

The only term that is potentially non-diagonal, ^(21) ’ therefore provides a good indication of how well- 
balanced the Gramians are, as we expect it to tend to the identity matrix. Similarly, we could investigate 
how close is to the identity matrix to investigate how well-balanced the controllability matrix 

is. 

We therefore choose to compute ||/ — V(2i)||| f, where the transformations have been computed 

using a set of snapshots corresponding to too = G = 20, 40, and 80 time units and applied to Gramians, 
which were computed for values of too = t2 ranging from 0 to 120. We only consider the top left (4 x 4), 


15 









Figure 5: (Color Online) Convergence of the balanced mode shapes, shown by plotting 

1 — \Ti(tyTi{t — At)I with respect to the final simulation time too used to compute T, for the first 5 
normalized balanced modes Tj = Here At = 0.2 and unstable modes are shown with thick 

lines. 



Figure 6: (Color Online) Balancing error of transformed Gramians computed from a simulation with final 
simulation time too = O and balanced using transformation matrices calculated from simulations with 
too = ti- The error is quantified by the top left (4 x 4) (thick), (8 x 8) (thin), and (12 x 12) (dashed) 
elements of ||/ — |C’(ii)C(2i)II |f for too = ti = 20, 40, and 80. The value of ti from each curve can be 
identified by the zero error values at t2 = ti. 

(8x8), and (12x12) elements of V(2i), to evaluate how balanced the resulting 4th, 8th, and 12th-order 
ROMs would be. 

Focusing first on the influence of ti, Fig. shows that in general, computing the matrix with a longer 
ti improves the balancing for all three orders of the ROM (for all except highest values of t 2 and close 
to ti = t2), as expected: theoretically, the higher the value of ti, the more converged the transformation 
matrices. Specifically, the three curves corresponding to ti = 20 have a large error regardless of the ROM 
order, as the transformation matrices are not converged. On the other hand, if ti is too large, the set 
of snapshots used to define the transformation does not allow an accurate identification of dynamically 
less significant modes. Therefore Gramians (and ROMs) which include these modes will not be properly 
balanced regardless of t2- In Fig.Eltherefore, the curve corresponding to a 12-mode Gramian with ti = 80 
always has a large balancing error. 

Gonsidering now the effect of t2, the quality of the balancing increases as ^2 approaches ti. At O = G 
the Gramians are exactly balanced by definition as mentioned above. For O > the error settles to a 
constant value, until eventually, a sharp degradation appears for the highest values of t 2 , for instance at 
t 2 ~ 100 for the ti = 80, 8-mode ROM. This behavior can be explained by the fact that the accuracy of 
the Gramians, AX'! and ZZ\ as opposed to the transformations, decreases if O is too large, as they are 
also approximated using state snapshots. 

In this section, we have shown that in practice, we cannot allow the system to fully converge by 
letting too +00, as this results in the loss of crucial information about the dynamically significant 
stable behavior of the system. However, the balanced modes do converge up to a critical too value and 
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Figure 7: (Color Online) ROM error as a function of the final simulation time too- Solid lines: snapshot 
projection method mm- Lines with circles: projection-free method. For both methods, the different 
lines correspond to ROMs with 2, 4, 6, 8, 10 and 12 states, from top to bottom. 


result in approximately balanced transformed Gramians. We therefore now turn our attention to the 
physical problem of the performance of reduced-order models obtained with this method, and compare 
them with ROMs obtained with existing methods. 

4.2.2 Analysis of the reduced-order models 

In this article, the Coo norm of the difference between the full-order transfer function and the reduced- 
order transfer function, defined as sup^gj^yg^ \G{iuj) — Gr{iuj)\ for SISO systems is used as a measure of 
ROM quality. In practice, we approximate this by: max,^g[_,^^ \G{iuj) — Gr{iuj)\ forwoo —t + 00 . Note 
that this transfer function norm is different both from the T-Loo error defined as sup^.j^gj-g^^o “ Gr-(s)|, 
which is infinite for unstable systems, and from the time-domain £ 00 (ft) norm usually used for signals 
x{t) and defined over some interval It as ess supfgjja:(t)|. The Coo error norm normalized by Coo 
norm of the full-order system will be referred to as the “ROM error” in the following paragraphs, i.e. 
||G — Gr.||oo/||G||oo- Potentially unstable systems can alternatively be compared using the (5i,-gap metric 
m, which can be interpreted as the “distance” between two plants from the point of view of their behavior 
in a feedback setting. However for the purposes of the present analysis, the Coo norm is an adequate 
measure of the difference between the two systems (in this case it represents the maximum distance 
between the loci of G{ioj) and Gr{iu})), as our focus is on model reduction as opposed to controller design 
in this article. 

Given the restrictions on the final simulation time mentioned in Sec. 13.21 and Sec. 14.2.11 we first 
investigate the quality of ROMs obtained with different too in Fig. [71 The behavior described in Sec. 
is apparent: as the simulation time increases the ROM error initially decreases. The lower bound on too 
related to the decay of the stable modes is illustrated with the black lines with no markers, corresponding 
to ROMs of different orders, but computed with the projection method of mm- The second lower 
bound on too, corresponding to the convergence of the unstable modes as they begin to dominate the 
impulse response appears to be slower, as shown by the lines with markers, corresponding to projection- 
free ROMs. Now, if too is sufficiently large for the less dynamically significant information to be lost, 
increasing the order of the ROM does not reduce the error any further. Additionally, each ROM order 
has a different optimal final simulation time: the error is minimized just before the unstable modes start 
to dominate the response enough to sharply increase the error. In the following paragraphs, an estimated 
value topt{r) for the optimal too corresponding to each ROM order r is used. These values are shown 
in Fig. El where it can be observed for this particular case that topt{r) seems to be decreasing roughly 
linearly as the model order is increased for r > 2. 

The existence of an optimal final simulation time can be further illustrated by plotting the singular 
values of the Hankel matrix (Fig. (TUI) for different values of too- Here, the HSVs corresponding to 
4th, 8th and 12th-order ROMs, obtained with the exact and snapshot projection methods as well as 
the projection-free method (using estimated optimal too values) are plotted. As too is increased, in the 
projection-free case, the HSVs corresponding to unstable modes grow, while the rest of the singular values 
tend towards their respective converged values. This explains the increasing precision of the ROMs as the 
hnal simulation time is increased, like in stable systems. On the other hand, the flat region of each curve, 
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Figure 8: (Color Online) Estimated optimal final simulation time topt{r) for different ROM orders r 
(symbols) and linear tend line plotted for r > 2. 
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Figure 9: (Color Online) ROM error as a function of the ROM order, obtained with too = topt{r) for the 
projection-free method. Crosses: exact projection method [50]. Filled circles: snapshot-based projection 
method mm- Open circles: projection-free method. 


where the singular values effectively have no dynamical significance, incorporates an increasing number 
of states as too grows. This again corresponds to the less significant information disappearing due to the 
unstable modes becoming too large, and can be used as a way to choose the appropriate ROM order, 
given an HSV distribution. For instance, with too = topi(4) = 108 (top line with diamonds in Fig. ITU)) , 
roughly 4 singular values are not in the flat region of the curve. As shown in Fig. [71 the error will start 
growing if too is increased further with r = 4. Conversely, increasing the ROM order to anything higher 
than 4 will not reduce the ROM error significantly when too — topi (4). 

In Fig. ini the error obtained with topt{r) is plotted as a function of the ROM order, and compared to 
the performance of the snapshot mm and exact jTD] projection methods. Clearly, the two projection 
approaches yield ROMs of the same quality. On the other hand, the projection-free method ROM error 
is of comparable order of magnitude, with the difference in the error increasing slightly with the ROM 
order. Note however that this difference is not significant: the error of a 12-state projection-free ROM is 
smaller than that of an 11-state ROM computed with the projection methods. 

For control design purposes, it is crucial that the ROMs are able to reproduce the full system’s impulse 
response and transfer function. In Fig. 1111 the impulse responses (a) and the transfer function gains (b) 
of the three methods are compared for different ROM orders. The three methods yield similar results for 
the three ROM orders considered. As the ROM order in increased in all three cases, the initial transient 
of the impulse response is estimated with increasing precision. Similarly, the high-frequency gain is better 
modeled with high ROM orders in the transfer function. 
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Figure 10: (Color Online) Hankel singular values of ROMs obtained with the exact projection method 
PP] (crosses), the snapshot-based projection method [TTl [T] (filled circles), and the projection-free method 
(open symbols) obtained with toe, = topt{r) for ROM orders of 4 (diamonds), 8 (squares), and 12 (open 
circles), from top to bottom. Note that the unstable singular values are not plotted for the first two 
methods as they are infinite by definition. 



Figure 11: (Color Online) Impulse response (a) and transfer function gain (b) of the full system (thick 
green dashed line) and the ROMs obtained with the exact projection method [PD] (thick blue dashed 
line), the snapshot-based projection method [Hill] (thin black dash-dotted line), and the projection-free 
method (thin red solid line). ROM orders: 4, 8, and 12, top to bottom for both the impulse responses 
and the transfer functions. 
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5 Application to a two-dimensional unstable flow-field 

The previous sections in this article focused on demonstrating the basis for using the method from a 
theoretical point of view and by comparing its performance to that of existing methods, when applied to 
the linearized complex Ginzburg-Landau equation. In this section, the applicability of the method in a 
more realistic scenario, i.e. for large-scale unstable flow-fields is investigated. The chosen test case is the 
flow over a cylinder at Reynolds number Re = UooD/v = 100, where Uoo is the incoming flow velocity, 
D is the cylinder diameter, and v is the kinematic viscosity. 

The forward direct numerical simulations run here are based on the immersed boundary projection 
method code developed by Taira and Colonius Si HD. The two-dimensional, viscous, incompressible 
Navier-Stokes equations expressed in vorticity form are discretized on a second-order accurate finite- 
volume staggered grid. A multi-grid algorithm is used, where the problem is solved on a series of nested 
identical Cartesian uniform grids, each of twice the physical extent and half the resolution of the previous 
one. A second-order Adams-Bashforth scheme is used to discretize the advection term(s) and a second- 
order Crank-Nicolson scheme is used for the other (linear) terms. For more information, the reader is 
referred to Si HI]. The linearized equations are therefore a discretization of: 

{ w = V X (uo X w) -I- V X (n X cuo) 

-Re-^Vx(Vxuj) + Vxf!B +Vxf, (20) 

ub= 0 , 

where w refers to vorticity, u to velocity, o to base flow quantities, fiB to immersed boundary forcing 
that imposes the no-slip condition on the cylinder surface, i.e. the second equation in Eq. (sni), where 
ub is the velocity on the body surface. The spatially distributed forcing function / is used to model the 
input (more details below). 

The corresponding continuous adjoint equations can be shown to be: 

{ —Cj^ = V X (wq X M+) — («■*■ X Uo) 

-i?e-iVx(Vxa;+)-b Vx//b-I-V x/+, (21) 

4 = 0 , 

where refers to adjoint quantities and /'*" forces the adjoint simulation at the sensor location (output 
location of the forward simulation), as opposed to /, which forces the forward simulation at the actuator 
(input) location. In this code however, the adjoint equations were obtained from the spatially discrete but 
temporally continuous linearized equations, where the time-marching scheme is self-adjoint in the present 
case since the base flow is the steady unstable equilibrium of the problem. Note that the same nested-grid 
method was used in both the forward and adjoint problem and this procedure is not self-adjoint so the 
adjoint equations we solve are technically not the exact discrete adjoint of the forward problem. The 
solver used for Eq. (I^Oll and (I21II is similar for both sets of equations, except for the advection terms and 
the fact that Eq. m runs backwards in time. 

The same grid as the one in |50j was used here, where 3 nested grids were used, and where the smallest 
grid has dimensions [—15,15] x [—5, 5] and the widest one [—60, 60] x [—20, 20]. Each grid has 1500 x 500 
square cells of length and width 6 = 0.02. The chosen input-output setup is also almost identical to the 
one in m- the flow is forced using a disk-shaped vertical body force of radius 1 (white disk in Fig. IT^. 
centered at (x, y) = (2,0) in the wake. The sensor measures the vertical velocity, also at (x, y) = (2,0) 
(black triangle in Fig. IT^ . The base flow, shown in Fig. [JU was computed using Selective Frequency 
Damping (SFD) [U [26j. The SFD parameters values chosen here were x = 0.4391 and A = 3.1974, as 
suggested by [2S| for the cylinder flow at Re = 100. 

In order to compute the reduced-order model, both the forward and the adjoint impulse responses 
were run for a long time (300 time units). Snapshots were stored every 0.2 time units and a Boole 
quadrature rule was used to scale the snapshots. Once the snapshots are stored, one can compute the 
Hankel singular values of the ROMs for different values of too, as shown in Fig.[T3| The general behavior 
is similar to the one in Sec. HJ but in this case the HSV distribution is more complex than in Fig. [ini 
Furthermore, both the transients and the long-term response of the system are of interest, it is less clear 
how one should quantify the ROM error. Nevertheless, since the cost of computing the SVD of the Hankel 
matrix is relatively small, one can find a compromise between ROM order and quality by trial and error 
by comparing the model impulse response to the one from the sensor. 

In this case, using only 125 snapshots, corresponding to a final simulation time of 25 time units, a 
12th-order model was obtained, whose impulse response is compared to the full system’s in Fig. [T3| It is 
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Figure 12: (Color Online) Unstable base flow over a cylinder at Re = 100, showing streamlines and 
contours of vorticity. The white circle represents the disk of diameter 1 where the input (a vertical body 
force) is applied and the black triangle shows the location of the output, a sensor measuring the vertical 
velocity at (2,0). 



Figure 13: (Color Online) Hankel singular values of ROMs obtained with the projection-free method with 
too = 13 (diamonds), too = 25 (squares), and too = 50 (open circles), from top to bottom. 


clear that the method successfully reproduces the impulse response in this case, both for times that are 
much longer than the final simulation time used to compute the model and for the initial transients. 

Note that the instability here is due to an unstable complex-conjugate pair of eigenvalues, i.e. the 
special case considered from a theoretical point of view in Appendix IA.3.21 This therefore demonstrates 
that the technique is also readily applicable in practice for this ubiquitous scenario. The HSVs corre¬ 
sponding to the two unstable complex conjugate modes clearly appear in Fig. [13] as the first two HSVs, 
which grow to large values as too is increased, and are always of the same order of magnitude. 

Another important point is that at t = 25, the output is |?/(t)|~ 1, which is only about 4 times larger 
than the first transient peak. In other words, we obtained a successful ROM of this unstable linear system 
without requiring an excessively long final simulation time or strong domination of the unstable modes. 
This suggests that with a sufficiently small impulse at the input, the method might be directly applicable 
to the initial linear response of some nonlinear systems. 


6 Conclusions 

In this article, we have shown that applying the snapshot-based, projection-free, approximate balanced 
truncation method to unstable systems theoretically yields a converged transformation that balances the 
Gramians for all too —^ +oo, including the unstable subspace. For stable systems, this method has been 
popular in recent years as it is straightforward to implement and scales well for large systems. The 
fact that it can be used without any modifications for unstable systems will allow the computation of 
BPOD-based ROMs and controllers of even three-dimensional unstable flow-fields, without the need to 
identify or project-out any global modes. 

Benchmarking the method using a one-dimensional, spatially developing, unstable system showed 
that it is not required for the system to converge to its theoretical limiting state for this approach to 
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Figure 14: (Color Online) Impulse response of the full cylinder flow (thick dashed line), compared to the 
12th-order ROM obtained with the projection-free method (thin solid line). The agreement between the 
two impulse responses is excellent, both for the initial transients shown in (a) and long term unstable 
response shown in (b). 


yield ROMs of similar quality to the ones obtained with existing methods [million]- In fact this limiting 
behavior cannot be reached exactly in practice due to finite-precision arithmetic: if the final simulation 
time becomes excessively large, only the unbounded growth of the most unstable modes will be conserved 
and the initial transient information from the impulse response will be lost. It was found that there is 
an optimal final simulation time value, which is a function of the desired ROM order, and methods to 
estimate this optimal value were discussed. 

In cases where the unmodified method yields unacceptable ROMs, or where the balanced modes must 
be identified more accurately, a method was outlined to circumvent the final time restriction. This method 
has a slightly increased cost as it requires projecting out some of the system’s unstable balanced modes, 
but their identification still does not require an Arnoldi solver. Additionally, if the system dimension 
is so large that storing a large number of snapshots is challenging, it was suggested that a (piecewise) 
Gaussian quadrature method (with varying sampling intervals) could be used when storing the snapshots, 
as opposed to a traditional Newton-Cotes approach with a constant sampling rate, as this is a simple 
way to significantly reduce the number of snapshots required to obtain a ROM of a given accuracy. 

Finally, the method was applied to a more realistic large-scale system: the supercritical (linearized) 
flow over a circular cylinder at Re = 100, and an excellent match was obtained with a I2th-order model. 
Only a short simulation and a small number of snapshots were required to obtain the ROM, where the 
final output amplitude was similar to the initial transients, suggesting that it may be possible to apply 
the method to the initial linear growth of nonlinear unstable flows. 

The standard approximate balanced truncation algorithm was therefore shown to be applicable to 
unstable systems both from a theoretical and a practical point of view. The implications of this are 
twofold. First, the simple numerical setup required to compute the snapshot-based balanced truncation 
of stable systems can be safely applied without any modifications to unstable systems. Second, the 
computation of balanced reduced-order models of unstable systems that are so large that other existing 
methods are not tractable is made possible with this approach. 

7 Acknowledgments 

This research is funded by the EPSRC. The authors would also like to thank Shervin Bagheri for allowing 
us to use the Ginzburg-Landau Equation MATLAB code, on which Sec.|3]is based, as well as Tim Colonius 
for the Immersed-Boundary Eractional Step code used in Sec. [S] 


22 






A Projection-free approximate balanced truncation of unstable 
systems: special cases 

In Sec. [21 it was assumed that the growth rates of the different eigenmodes of the flow field were dis¬ 
tinct, so that one eigenmode would eventually dominate the impulse response. The purpose of this 
section is to show that the proof extends in a similar way when some of the eigenmodes have the same 
growth rate. The general case where the nth to (n -I- m)th eigenvalues have the same growth rate: 
Re(A„) = Re(A„+i) = ... = Re(A„+m) = a, but the imaginary parts Im(A„) = oJn differ: ^ uin+i ^ 

... ^ U!n+m is considered first and the special cases of repeated eigenvalues and complex conjugate 
eigenvalue pairs are then treated separately. 


A.l General case 


A.1.1 Controllability and observability singular vectors 

Let us assume that the (n — 1) most unstable modes are projected out as in Sec. [S] Let us also write 
n = diag([ (jJn ... oon+m ])■ The impulse response state x{t) = becomes: 


lim Tcn-lx{t)=Tcn-l\pn 


lim Tqyi—iX — Tq^iti—i \Pn 

too—^ + 00 ^ 


= [ 


'^cn 


^cn-\-m 


^ cn 


iVtt 


Pn-\-m \ ^ 


• • • Pn-\-r 


Be°‘\ 






^ cn-\-m 


- 



cn 


t 

_ 

. ^cn+m _ 


The left and right singular vectors of Tcn-iX thus tend to the subspaces spanned by Tcn-i [Pn ■ ■ ■ Pn+m ] 
and [ /3n ... Pn+m ] respectively. Ton-iZ behaves in an analogous way and hence its left and right 

singular vectors tend to the subspaces spanned by Ton-i [ Qn ■ ■ ■ Qn+m ] and [ G ■ • ■ ^n+m ] re¬ 
spectively. However [ Ucn ■ ■ ■ Ucn+m ] and [ Uon ■ ■ ■ Uon+m ] do not generally converge to con¬ 
stants due to the oscillatory term . 


A.1.2 Hankel matrix 

Next, as in Sec. O we consider the projected Hankel matrix: 


lim Tin-iMuTm-i = lim Tin-iZ"^ XTm-i 

too—>• + 00 too—> + 00 


= [ 


Xin • • • ^n+m 





1 - 

1 _ 


.yt 

_ ^n+m _ 


— [ '^on ■ ■ ■ '^on+m ] ^ 

where: 


“^cn+m 


M = 


f^on+m'^on+m _ 


'^cn^cn • • • ^cn+m^cn+m • 


( 22 ) 


The left and right Hankel singular vectors therefore tend to the subspaces spanned by [ Von ■ ■ ■ Xon+m ] 
and [ Vcn ■ ■ ■ Xcn+m ] respectively, which are already (individually) orthonormal bases by definition. 
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We thus only need unitary matrices to identify the Hankel singular yectors and values, which we can get 
from the SVD of M: 


M 


Ro 



Therefore in general: 

^ Un ■ • ■ Un+m ] “ [ ^on 
^ Vn • ■ • Vn+m ] “ [ '^cn 

where Rc and Ro are not necessarily constant. 


Rl 




'^on+m ] R'Oi 
'^cn-\-m ] Rc: 


(23) 


(24) 


A. 1.3 Balancing transformations 

For large too, only the nth to (n + m)th columns of the transformation matrix T = must 

be considered since Eq. (l24l) implies that vl^Vj = 0 ii n < i < n + m and j ^ [n,n + m], and hence the 
rest of T is independent of these (m + 1) columns. Using Eq. (IMl) . these (m + 1) columns become: 


lim [ r„ ... Tn+m ] 

too —¥-\-00 ^ 




^cn 



r - 1/2 

<Jn 


[ '^cn 



^ cn-\-m 

Rc 


- 1/2 

^n+m 


Therefore, in general, the subspace spanned by the balanced modes Ti for n < i < n + m converges to 
the subspace spanned by Tcn-i \ Pn ■ ■ ■ Pn+m ] ■ Similarly, it can be shown that the adjoint balanced 
modes Si for n < i < n + m must span the same subspace as Ton-i \ Qn • • • <ln+m ]. 


A. 1.4 Transformed Gramians 


Finally, since Eq. m is still valid here. Wo is diagonalized by T if V( 2 i) = We again only need to 

consider the nth to (n + m)th columns of each matrix since the other columns are unaffected by these 
modes for large ti and ^ 2 - Using Eq. (El): 

V-o, . . . '^n+m ] ^21) ~ I- ' * ' ^cn+m ] {^l')Rc{2l): 

^ Vn . . . Vn-\-m ] — [ ^cn ■ • ■ "i^cn+m ] Rc{ll): 

=^[ '^n Vn+m ] ( 2 I) [ Vn+m ] -^c(21)-^c(ll) • 

In general, therefore becomes: 


^(n)^( 2 i) 


'I 0 O' 

0 ^I(ll)-^c( 21 ) 0 

0 0/ 


The transformed Gramian T{ti)'^Wo{t 2 )T{ti) is thus fully diagonal, except for (m + 1) columns which 
have a dense ((m + 1) x (m + 1)) block along the diagonal. An analogous derivation applied to the 
controllability Gramian leads to the conclusion that the Gramians are balanced (diagonal and equal), 
except for the (m + 1) columns corresponding to the modes with equal growth rates. The fact that the 
transformation does not in general balance unstable modes with identical growth rates is not problematic 
because, as noted above, the subspace spanned by the (m + 1) balanced modes does converge, and the 
rest of the system remains unaffected. In other words, the (m + 1) modes are balanced as a whole. As 
it is undesirable to truncate unstable modes for control purposes, the transformations T and S are still 
adequate transformations as long as they are computed with a set a snapshots that allows this subspace 
to converge. 

On the other hand, if Rc and Ro were to tend to constant matrices, the Gramians would become 
fully balanced by T and S for large too- We will focus on special cases where this happens in the sections 
below. 
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A.2 Repeated eigenvalues 

In this case, A„ = A„+i = ... = Xn+m = a + iiu. In order to analyse the singular values and vectors of X 
in this case the transformation is applied to the full controllability Gramian tbc(too) = /q” xx^dt = XX^-. 

Tcn-lWc{t^)Tcn-l = f Tcn-lXX^Tcn-ldt, 

Jo 


which gives: 


^+oo ^cn—l^^c(^oo)^cn—1 — [ '^cn ■ ■ • '^cn+m ] 


'o-cn 


r 7 -t j 

^cn 

2 


t 

^cn+m . 


_ ^cn+m _ 







pi 

Tcn—1 \_Pn ■ 

• Pn-\-m ] 


BB\qn . 

• Qn-\-m ] 




_ ^n-\-m _ 



_ Pn-\-m _ 


T ,„_1 


(25) 


since A + A* = 2a. The term to the left of /g“ e‘^°'*dt in Eq. (1^51) is just a constant matrix and the 
integral itself is a scalar. As a result [ Ucn ■ ■ ■ Ucn+m ] converge for large too and e'^°‘*dt = 
(e2“‘“ ~ l) / (2a) ~ e2“*“/ (2q;). The singular values can thus be written ~ d-cie°'*°° I where 
are the singular values of the constant part of Eq. (l25ll . If a = 0, then the modes are marginally stable 
and dt = too so a^i ~ Analogous conclusions can be drawn regarding the observability 

singular values and vectors: [ Uon ■ ■ ■ Uon+m ] converge and tJoi ~ S-oie°‘*°° j (and aoi ~ doi\ft^ 
for a = 0). 

The SVD of the matrix M defined in Eq. fiTIh can thus be written: 


where 


lim M 

too —t-\-00 


Ro 





2a 


Ro 








M 


(Jon^ 


t 

on 


_ ^on-\-m'^on-\-m 


'^cn^cn 


^cn+m^cn+m 


(26) 


Here ai = f {2a) (and Ui = ditoo for a = 0), and M tends to a constant matrix (and hence so do 

its singular values (t„, and vectors Ro, Re)- 

For large too, the nth to (n + m)th columns of the transformation matrix become: 





'--1/2 


Rc 


cn+m 


~-l /2 


^n+m 


Since Rc and Ro converge to constants for large too, we can conclude that for repeated eigenvalues, the 
matrices T and S converge and fully balance the Gramians. 
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A.3 Complex conjugate eigenvalues 

If again there are (n— 1) modes that are more unstable than the complex conjugate pair, then we proceed 
to identify all of them and project them out, as in Sec. [31 


— '^cn—1 \_Pn ■ • ■ Pua; ] 
lim TQYi—\x{t'j — 1 [ Pn 

t^ + OC 


o^n.i 




7*t 




■gA„t 

0 

\ 1 

o 

_1 

eKt 



B. 


Finally, this results in; 


lim Tcn-lWc{too)Tcn-l=l 

too^ + OO / 


Tcn—l \_Pn 




Pn'^cn—l 


q*-^BB'<q*e^°‘* _ 

_ Pn^Tcn-l_ 


A.3.1 Marginally stable case 

In this case = 1, and the projected Gramian simplifies to: 


where: 


lim Tcn-lWc{too)Tcn-l 

■ oo—^-\-00 



q^BBU*e+' 

Pn'^cn—l 

_q*'^BB^qe- 


Pn^Tcn-l_ 


^ ^{±2iuJntao) _ I 

dz2zCc^77,toO 


lim = 0, 

too—¥-\-00 


and hence the projected Gramian can be written: 


lim Tcn-lWc{tac)Tcn-l 

t—)- + oo 

Pn^cn —1 

P*^Tcn-l 

= tooWc, 


— f T 

— brxy J- r.n.— 


oo-^ cn — l\ Pn Pn 


q^BB^ 0 
0 q*'!BBU* 


where Wc is constant, and hence has constant singular vectors [ Ucn Ucn+i ] and singular values that 
can be written: toc^cn and toodcn+i, with den, ^cn+i constant real scalars. The observability Gramian 
behaves analogously. As a consequence, the matrix M defined in Eq. (12^ can be written M = tcoM, 
where M is constant, so Rc and Ro tend to constants and Ui = tao^i, where di are the constant singular 
values of M, just like in Sec. IA.2I The conclusions that T and S tend to constant matrices that fully 
balance the Gramians follow in the same way. 


A.3.2 Unstable case 

If a > 0, we obtain: 


lim Tc„_ilTc(too)T'c„_i 

t—^ + oo 


Tcn-l[Pn P*n] 


Unlike the marginally stable case however: 
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Therefore, even for large foo the left controllability (observability) singular vectors oscillate but stay in 
the two-dimensional plane defined by the eigenvectors. Similarly the Hankel singular values grow on 
average exponentially as , but the growth rate also oscillates around this mean trend. As discussed 
in Sec. EH for each unstable complex conjugate pair of modes, two columns of T and two rows of S 
will not converge regardless of too (but will stay bounded), and a full (2 x 2) block will be present along 
the diagonal of the transformed Gramians, while the rest of the matrix will be independently balanced. 
If too is large enough for the subspace spanned by the two corresponding columns of T (and rows of 
S) to be converged, then any T and S transformations can therefore be used an adequate balancing 
transformations, where the complex conjugate pair is considered to be balanced as a whole. 

B Induction proof for the Hankel matrix SVD 

As in Proposition 1 in Sec.[31 we will use a proof by induction to show that for large too, the singular vectors 
and values corresponding to unstable modes of the Hankel Matrix Mh — Z^X = = Vo'SoU^Uc'Z'cV^ 

tend to: 


{ ^i — '^oi, 

Vi = Vci- 

We choose the induction hypothesis (I„) to be the fact that the proof holds for alH < n. For the base 
case (Ii), if the first mode is unstable: 

lim Mh = VoiUoiul^UdCfcivli = uicriuj, 

too—^ + <30 


since, if too is large, CTci Ud and cToi Uoi for i > 1. (Ii) therefore holds: ui = Voi, cri = (ToluolZero’d 
and vi = Vci- 

In the inductive step, we prove that if (I„) holds for some rank n, then (I„+i) also holds. The 
transformation matrices Tu defined in Eq. (nm project out all the left singular vectors Uj for all j < i, 
i.e. TiiUj = 0, TijUi = Ui. acts in an analogous way on the right singular vectors. As we are assuming 
(In) to hold for rank n: 


IlnXlnlrn — ['^n-t-1 ] 


= Tir, 


^n+1 


'^n+1 


^on+1 
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O’onUon-l-l 
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^cn+1 

_ _ 


. O’orix ^orix - 


. CTcrix 'aJrix - 
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As a result, if the (n -|- l)th mode is unstable: 

lim TinMnTrn 

too ^ + CX3 

because if t^o is large, ad ^ acj and aoi ^ aoj for i < j. Since Vo and Vc are individually orthonormal 
bases, Von+i is normal to Voi = Ui i < n and fcn+i is normal to Vd = Vi i < n: 

lim TinMHTrn = '^on+imon+ltto.^_l_X'Ucn+lf^cn+l) “^cn+l ’ 

too—>- + 00 V / 

and therefore (T„+i) holds too: u„+i = Von+i, o-„+i = (Jon+iv^o^^iUcn+iUcn+i and Vn+i = Vcn+i- This 
completes the inductive step and, along with the base step, this concludes the proof by induction. 
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C Proof of the convergence of the ratio of singular values 


In this section, we prove that o'ci/®'* tends to a constant for large too- 
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Only ||/3j|| and ||^i|l are not constants for large too in the above expression and: 
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We can therefore use the approximation: 
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where we have defined Ai + A* = 2ai and ai G M. Similarly: 
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